Analysis of microtubule motion due to drag from kinesin walkers 
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We analyze the nonlinear waves that propagate on a microtubule that is tethered at its minus end 
due to kinesin walking on it, as is seen during the fluid mixing caused by cytoplasmic streaming in 
Drosophila oocytes. The model we use assumes that the microtubule can be modeled as an elastic 
string in a viscous medium. The effect of the kinesin is to apply a force tangential to the microtubule 
and we also consider the addition of a uniform cytoplasmic velocity field. We show that travelling 
wave solutions exist and analyze their properties. There exist scale invariant families of solutions 
and solutions can exist that are flat or helical. The relationship between the period and wavelength 
is obtained by both analytic and numerical means. Numerical implementation of the equation of 
motion verifies our analytical predictions. 
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I. MOTION OF MICROTUBULE WITH KINESIN WALKERS 

Here we study the motion of microtubules in the context of cytoplasmic streaming in stage 10B-11 Drosophila 
oocytes. At that stage, the roughly hemispherical oocyte is bounded by a plasma membrane and an underlying cortex 
comprised of an actin filament meshwork. Long microtubules, with minus ends attached to the cortex, have their 
plus ends free. Kinesin-1 motor proteins that walk along the microtubules toward plus ends generate opposing forces 
on the microtubule and the surrounding viscous cytoplasm. Cytoplasm is thus moved toward plus ends, generating 
vigorous flows for mixing. Because minus ends are anchored, the equal but opposite force on a microtubule generates 
dynamic bending along its length. The results of this work have been used to understand this phenomenon in recent 
work by the authors 1 1 . 

We begin with the equation for a microtubule with an applied force tangent to the direction of chain. The position 
of the microtubule at arclength s and at time t is denoted r(s,t). As is usual at small scales, all inertial effects are 
negligible and the system is dominated by the drag coefficient per unit length i/, 



at as 4 os as as 



C denotes the elastic bending constant of the microtubule. The tension T(s) enforces the inextensibility of micro- 
tubules which can be written as |<9r/<9s| = 1. Kinesin motors walking along the microtubule are assumed to be 
carrying cargo that acts as an impeller pushing the surrounding cytoplasmic medium. By Newton's third law, this 
force produces a drag on the impellers that transmits this force to the attached microtubule. Assuming a high density 
of kinesin walkers, the effective force on the microtubule will be in the direction of the local tangent to r(s) which is 
dr/ds. The coefficient gives the strength of the force. Lastly we add a uniform flow field representing the velocity 
of the cytoplasm v s in the vicinity of the microtubule which we take to be in the k direction. 
We rescale to dimensionless variables 

cr = s/p , u = r/po, t = tuj , T' = Tpl/C, and h = vv s /f k (2) 



so that 



du d 4 u d , du du - 
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This requires that 



w = fk/(vpo), and pi = C/f k . (4) 



Po and wo are constants that make our new variables dimcnsionlcss. 

We will consider solutions where one end (the minus end) is tethered to a point, say the origin, so that r(s = 0, t) = 
and is freely hinged at that point. The other end is free. The total arclength of the microtubule is denoted by L. 
The general characteristic of all numerical solutions to Eq. Q] is that for large s the solutions appear to be travelling 
waves. In fact, the form of solution becomes independent of L so in fact we could consider this problem in the limit 
L — > oo. For example we find solutions that asymptotically become helical with a fixed path and radius for large s. 
Other solutions are planar and are periodic in s. In all these cases the solution has the form of a travelling wave that 
we analyze in detail below. 

II. STEADY STATE TRAVELLING WAVE SOLUTIONS 

Now we look for steady state travelling wave solutions. First we clarify what this means. The whole microtubule 
should not be translating because one end, (far away from where the travelling wave solution is valid) is tethered. 
The position averaged over one period of a monomer, (u), should be time independent. If the microtubule is stretched 
out in one particular direction, there should be a displacement of the microtubule about a straight line solution, 

u(cr, t) = \id{<7 — vt) + aak. (5) 

The term describes a displacement that travels along the backbone of the chain maintaining its shape. Therefore 
we require 

d(u d ) 



(6) 
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Substituting Eq. [5] into Eq. Ogives a left hand side of —vdu.d/da and a term of the form —d\id/da on the right 
hand side. Choosing v = 1 and in addition a = h, we are left with 

We define w = du/da and note that the inextensibility requirement \dr/ds\ = 1 implies that |w| = 1. We can write 

w = — h ak (8) 

oa 

Note that w(cr, r) depends only on a — r. Because w is a function of only one variable £ = a — r, we can integrate 
Eq. [7] with respect to £ obtaining 

^--T'(a)w = -gk (9) 



The integration constant on the right hand side —gk must lie in the k direction by symmetry. This is the same 
equation as that of the classical mechanics of a particle of unit mass, with position w travelling on the surface of a 
unit sphere under the influence of gravity, with the variable £ being analogous to time. The term — T" is a normal 
force that constrains the particle to stay on the sphere's surface. 

From Eq. [6] and Eq. [8] we see that the motion is subject to the constraint 

(w) = ak (10) 

The constant g and the initial conditions of the particle must be chosen so as to satisfy this constraint. There is an 
infinite number of solutions satisfying these conditions leading to different steady state solutions for the microtubule. 

The spherical pendulum Eq. [9] has been analyzed [2[ in detail. In general the orbits will not be closed, but we will 
relegate our discussion below to two cases where this is the case, circular motion in Eq. [5J corresponding to helical 
conformations of a microtubule, and motion in the x — z plane, corresponding to solutions similar in shape to cycloids. 



A. Scale Invariant Family of Solutions 

One important point to notice about the general structure of Eq. [5] is its invariance under change of £, that is 
£ — > A£. A change in scale of A does not modify the solution, because T"(£) and g are both constraints, and can 
therefore be chosen arbitrarily. So that if w(£) is a solution to Eq. [9j so is w(A£). Furthermore all solutions of 
this form will have the same (w) and hence the same a. Note that a length scale shift of w = wo(A£) corresponds 
to a displacement of u = (l/A)uo(A£), which changes scale because of the 1/A prefactor in this expression. This 
means that for every steady state conformation of the microtubule, there exists a family of steady state solutions with 
identical shape but arbitrary size. 



B. Circular Orbits 

Circular orbits are solutions to Eq. [SJ Eq. [TU] requires that the vertical height of the orbit be at w z = a. Therefore 
the radius perpendicular to k (the w x , w y plane) is 

w± = s/l - a 2 . (11) 

It is first easiest to analyze this by direct analogy to the classical mechanical problem of a unit mass particle moving 
on a unit sphere (the spherical pendulum). If the height above the sphere is a, then application of Newton's laws 



gives that the tension must balance the force of gravity g and the centripetal acceleration v 2 /w± giving 

g _ a 



(12) 



v 2 1 — a 2 

Where v is not the real velocity of the microtubule but \dw±/d£\. This is the rate at which the tangent vector rotates 
in the x — y plane. Because \d\id/d£\ = |w^| = r c is a constant, w<j is a circle. Therefore according to Eq. O this will 
give a helical conformation of the microtubule. 

In the spherical pendulum analogy, the angular velocity of the particle is Q p = v/w±. The total arclcngth of the 
chain (in dimensionless units) of a single period of the chain is 

S = -^2= (13) 

and Sflp = 27r so that S = 27r/f2 p = wj_/v. Comparing this with Eq. [13] we obtain 



R= — \l\~a 2 . (14) 
v 

Using Eq. [IT] this gives 

\-a 2 l-a 2 

R = —, = (15) 

1^1 

Because g is a free parameter, v can be any positive number and therefore Eq. II 51 implies that the radius of the helix 
can be arbitrary. 

Because the velocity of this travelling wave is unity, the spatial and temporal periodicities are equal. Therefore the 
period of rotation of the helix P' is equal to S in Eq. I13[ in the dimensionless units that we are using. Therefore the 
relationship between the period and the pitch is 

2ttR , . 

p = -Jr=^i (16) 

V 1 — Or 

In terms of our original units the period is 

fkVl - a 2 
C. Two Dimensional Solutions 

We now consider orbits in the x — z plane. This corresponds to a pendulum swinging around vertically through the 
bottom of the sphere. This is a standard introductory mechanics problem. Letting 9 denote the angle with respect 
to the — k direction, 



(9= -9 sin A (18) 
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where the dots that normally represent a time derivative are really derivatives with respect to £. Low energy solutions 
where the total energy E is less than g give 9 oscillating between two bounds. In this case the curve for the 
microtubule will look sinusoidal. Beyond the separatrix where E > g, 9 will increase without bound. This gives rise 
to microtubule conformations that look close to prolate cycloids. In other words, the microtubule forms a travelling 
wave that periodically loops backwards. This is shown in Fig. Q] 




Z 



FIG. 1. A two dimensional travelling wave solution for a microtubule that looks close to a prolate cycloid but differs from it 
in functional form. 



Small values of \a\ correspond to low values of g/E. In this case, the particle spends almost equal times at all 
points on the circle giving a small value of (w) = a. As g/E increases to 1, a increases and the loops become tighter; 
a situation that is very costly energetically. However this is not the only solution to these equations for a fixed value 
of a. If instead we consider solutions with E < g, then large g/E corresponds to small oscillations of w about —k. In 
this case, |a| is close to 1 and the curve has no tight loops. The shape of the microtubule is close to a sine wave. 

To obtain the relationship between G = g/E and a we use conservation of energy to obtain the period T, 



9 = ^/2E(l + Gcos9) (19) 

The time average of w in the k direction is 

r cose — de 

(cos6>) = -pr^ j — (20) 

n / ; 0-9 

J0 y/l+Gcose) 

This can be expressed in terms of elliptic integrals as 

w = fi (I_(1+G) Wtl ) ( } 

Where £ and K, are the complete Elliptic integrals of the first and second kind respectively. For bounded orbits, below 



the separatrix, the corresponding expression can be calculated giving 



A*?) 



A plot of the a = (cos 9} versus 1/G = E/g is shown in Fig. [2] using these two expressions, a 
G — > 0. Note that there are two possible values of E/g corresponding to one value of positive a. 



(22) 



in the limit as 



a o.o 




FIG. 2. A plot of the value of a as a function of 1/G. 



The case of a — corresponds to the case of circular orbits discussed in Sec. Ill Bl The cytoplasmic streaming 
velocity is oc h (see Eq. [3]). When h — 0, we have shown that the solutions are circular. If h is now made very small, 
we expect that the solution will only be slightly perturbed from circular solutions. This corresponds to a small value 
of g/E. Therefore the value of G that is chosen for small enough h should correspond to the larger value of 1/G 
shown in Fig. [21 This corresponds to almost circular loops shifting slightly forward after every turn. As h increases 
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the loops become tighter giving rise to a high elastic bending energy. At some point therefore, we expect a transition 
to another state. In fact, simulations discussed in Sec. MI Al show that near a flat surface, the microtubule transitions 
out of the plane to a helical shape at h ps .385 for 128 links in a chain. For still larger h, the microtubule becomes 
flat again looking close to sinusoidal. 



III. NUMERICAL IMPLEMENTATION OF THE FULL SOLUTIONS 



By comparing the steady state solutions found above for large arclength s and time t, we will see that these solutions 
are physically realistic ones to consider. As we will see, the full equations of motions go to solutions of this type. 
However this is not a complete description of the problem, because of these solutions a particular one is selected 
from a whole family of solutions. The same behavior occurs in other problems in pattern selection such as dendritic 
growth [3j, |4|. 



A. Comparison With Simulation 

Eq. [I] was analyzed numerically using a method similar to that used earlier in the context of gel electrophoresis 5j,|6 
Link length drift was handled using a similar procedure to that implemented for chains with inertia [7[ One end was 
constrained to have coordinates at the origin, while the other was free. 

A Runge Kutta time step was 0.001, C = 20, = 2 and simulations were performed with N — 64 or TV = 128 
links, as will be noted below. 

We first consider the case of a microtubule with no wall or other external forces aside from the external velocity 
field. As a function of the rescaled external velocity field h, the radius of curvature and period were calculated once the 
system had reached steady state. Using rescaled variables as defined by Eqs. [5] and 0] we can write the dimensionlcss 
radius of curvature as R' c = R c / po, and the dimensionlcss period as P' = coP. The radius of curvature was calculated 
at the middle of the chain to reduce finite size effects. The results are shown in Figs. [3]and|U The data for R c , Fig. 
[3] using 64 links, is very close to those of 128 links except for the highest h values. The data for the period, Fig. QJa) 
show good agreement for both chain sizes at all values of h studied. 

We now test the analytic predictions that we made relating the period to the radius of the helix and a given by 
Eq. 1161 and our conclusion that a — h. The relationship between the radius of curvature R c and the radius R of the 
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FIG. 3. The rescaled radius of curvature R c / po versus the rescaled velocity field h = a. The green triangles are for 64 link 
chains, the red + symbols are for chains of 128 links. 

helix is readily calculated to be R = i? c (l — a 2 ). Therefore Eq. [TBI gives 

P' = 2nR' c Vl - a 2 (23) 

With this prediction and using the data for R c in Fig. [3J we can independently calculate P' for 64 link chains. We 
can only do this where finite size effects are not important and therefore omit the highest two h values from this 
analysis. The data, Fig. 0|b), show excellent agreement to within the differences expected by finite size effects. This 
corroborates our analytical analysis of steady state solutions. 

The case of h = a = displayed in Fig. [3] can be written in terms of the original dimensional variables of Eq. [1] 



10 



(a) 




(b) 



18 
17 
16 
15 
o. 14 
13 
12 
11 



0.4 0.6 
h 



FIG. 4. (a)The rescaled period P' = Pui versus the rescaled velocity field h — a. The blue circles are for 64 link chains, the 
red + symbols are for chains of 128 links, (b) The rescaled period P' — Poj versus the rescaled velocity field h = a calculated 
by using Eq. [23] and data on the radius of the helix versus h shown in Fig. [3] The result is shown by the green triangles. The 
blue circles are the same data shown in Fig. [4] by directly measuring the period. 



using Eq. [4] 



R c = (C/(/3/ fc )) 1/3 - 



(24) 



From the numerical solution, this gives (3 = 0.05 ± 0.0005. This is useful in analyzing experimental data, as R c is 
nearly constant for h < 0.7. 

Next we consider the presence of a wall. We introduced a force representing a wall in the y — z plane of the form 

(25) 



w {x + lf 

which is only present for x < 0. The force is singular at x = —1 preventing the chain from crossing that plane. The 
period as a function of h is shown in Fig. [5]for C = 20 and chains with 128 links. 

For small values of h, the period does not vary much. The shapes obtained are flat and appear to be the same as 
those found in III CI However there is non-analytic behavior at h sa 0.385 corresponding to a transition to flattened 
helical waves. Then at h « 0.425 the solution becomes almost completely flat and showing waves that look closer to 
sine waves, corresponding to the flat sinusoidal-like solutions studied in III CI 

We now consider the dependence of period and curvature on the length of chains L. Our analytical solutions have 
been in the limit that the chain length L — > oo. For finite length chains, we expect corrections to this behavior. This 
is important to analyze in relationship to experiments using microtubule gliding assays [8J. In this case we consider 
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h 

FIG. 5. The rescaled period P' — Puj versus the rescaled velocity field h — a for chains of 128 links. 



h = and study how the rotating spiral wave solutions vary with increasing length. Fig. [BJa) shows the dependence 
of the radius of curvature measured halfway along the arclength, as a function of chain length. As usual, the rescaled 
dimensionless variables, see Eqs. [5] and HI have been used. The curve shows a non-monotonic dependence on L that 
levels off at L as 20. Fig. [6tb) shows that the rescaled period is slightly non-monotonic as well but decreases to a 
constant value at L k 15. In comparison with gliding assays, we took [l| L = 14.0, where values of parameters are 
within 3% of their asymptotic values. 

Note that in Fig. [^b) the period diverges for finite L ss 3. Below this point, the solution is a static straight line. 
This is similar to the usual buckling transition for a finite length elastic rod. The microtubule needs to be sufficiently 
long to undergo the dynamical instability analyzed here. 
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FIG. 6. (a) The rescaled radius of curvature measured at the middle link of a chain versus the rescaled chain length in steady 
state for h — 0. (b) The rescaled period versus the rescaled chain length in steady state for h = 0. 



B. Selection of Scale 

As noted in section Hi A[ there exists a continuous family of steady state solutions because if w(£) is a solution, 
then so is w(A£) for arbitrary A. The general way that a particular value of A is selected is likely to be due to the 
same mechanism as in other pattern formation problems 3|, |J| . The travelling wave solutions ignore the boundary 
conditions at the ends s = and s = L. Far from those ends, we have found a continuous family of solutions. 
However these solutions become invalid close to the ends. The travelling wave solutions are only valid in the limit of 
<< s << L. The full solution must match to one of these travelling solutions in that region but will differ greatly 
near the ends. 

As an example, consider the circular solutions of Sec. Ill Bl with h = a = v s = 0. There the steady state solution is 
a rotating circle wrapped around on itself of arbitrary radius R. However near the end s = 0, the solutions goes into 
the circle's center. Similarly it deviates from a circle at s = L. In analogy with other pattern growth problems such 
as the "Geometric Model" 3|, where the mathematics have been analyzed in detail, we expect that the boundary 
conditions imposed on the ends, will only be satisfied for certain values of R. If there is more than one allowed value 
of R, the value picked out will be the most stable of these. In the Geometric Model the form of the equations is 
much simpler, making it possible to understand the overall structure of the problem more easily. In the present case, 
a precise understanding of the numerical results will be the subject of future research. 
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IV. CONCLUSIONS 

We have analyzed a model that describes a microtubule tethered at its minus end and subject to viscous drag and 
forces due to kincsin walking toward the plus end. The walkers provide a force tangent to the local direction of the 
microtubule and cause an instability that causes nonlinear waves to travel along its arclength. The scale of the wave 
can be arbitrary. An open question now is to find an analytic model that predicts the precise scale that is selected. 
There are many other questions that are interesting. What is the motion when there are obstacles that prevent free 
motion of the microtubule? Do these lead to chaotic behavior? What happens in the case of many microtubules 
interacting through steric and hydrodynamic interactions? This last question is difficult to answer and will be tied 
into the question of large scale collective motion of forests of microtubules and the velocity flows that they produce. 
Velocity flows giving rise to micro-mixing have been studied in related systems 9Ml2l] and it would be very interesting 
to understand the large scale cytoplasmic motion for Drosophila oocytes. 
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Note: After the completion of this work, we came across work of Bourdieu et al. [13| which analyzed their experi- 



mental results on spirals in mysosin and kinesin gliding assays using the two dimensional version of Eq. [T] obtaining 
the identical scaling, and simulating it using a nearly identical approach. This provides further evidence that the 
mechanism we have proposed for cytoplasmic streaming is physically viable. 
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